home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dtrevc.f < prev    next >
Text File  |  1996-07-19  |  35KB  |  990 lines
  1.       SUBROUTINE DTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
  2.      $                   LDVR, MM, M, WORK, INFO )
  3. *
  4. *  -- LAPACK routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     September 30, 1994
  8. *
  9. *     .. Scalar Arguments ..
  10.       CHARACTER          HOWMNY, SIDE
  11.       INTEGER            INFO, LDT, LDVL, LDVR, M, MM, N
  12. *     ..
  13. *     .. Array Arguments ..
  14.       LOGICAL            SELECT( * )
  15.       DOUBLE PRECISION   T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
  16.      $                   WORK( * )
  17. *     ..
  18. *
  19. *  Purpose
  20. *  =======
  21. *
  22. *  DTREVC computes some or all of the right and/or left eigenvectors of
  23. *  a real upper quasi-triangular matrix T.
  24. *
  25. *  The right eigenvector x and the left eigenvector y of T corresponding
  26. *  to an eigenvalue w are defined by:
  27. *
  28. *               T*x = w*x,     y'*T = w*y'
  29. *
  30. *  where y' denotes the conjugate transpose of the vector y.
  31. *
  32. *  If all eigenvectors are requested, the routine may either return the
  33. *  matrices X and/or Y of right or left eigenvectors of T, or the
  34. *  products Q*X and/or Q*Y, where Q is an input orthogonal
  35. *  matrix. If T was obtained from the real-Schur factorization of an
  36. *  original matrix A = Q*T*Q', then Q*X and Q*Y are the matrices of
  37. *  right or left eigenvectors of A.
  38. *
  39. *  T must be in Schur canonical form (as returned by DHSEQR), that is,
  40. *  block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
  41. *  2-by-2 diagonal block has its diagonal elements equal and its
  42. *  off-diagonal elements of opposite sign.  Corresponding to each 2-by-2
  43. *  diagonal block is a complex conjugate pair of eigenvalues and
  44. *  eigenvectors; only one eigenvector of the pair is computed, namely
  45. *  the one corresponding to the eigenvalue with positive imaginary part.
  46. *
  47. *
  48. *  Arguments
  49. *  =========
  50. *
  51. *  SIDE    (input) CHARACTER*1
  52. *          = 'R':  compute right eigenvectors only;
  53. *          = 'L':  compute left eigenvectors only;
  54. *          = 'B':  compute both right and left eigenvectors.
  55. *
  56. *  HOWMNY  (input) CHARACTER*1
  57. *          = 'A':  compute all right and/or left eigenvectors;
  58. *          = 'B':  compute all right and/or left eigenvectors,
  59. *                  and backtransform them using the input matrices
  60. *                  supplied in VR and/or VL;
  61. *          = 'S':  compute selected right and/or left eigenvectors,
  62. *                  specified by the logical array SELECT.
  63. *
  64. *  SELECT  (input/output) LOGICAL array, dimension (N)
  65. *          If HOWMNY = 'S', SELECT specifies the eigenvectors to be
  66. *          computed.
  67. *          If HOWMNY = 'A' or 'B', SELECT is not referenced.
  68. *          To select the real eigenvector corresponding to a real
  69. *          eigenvalue w(j), SELECT(j) must be set to .TRUE..  To select
  70. *          the complex eigenvector corresponding to a complex conjugate
  71. *          pair w(j) and w(j+1), either SELECT(j) or SELECT(j+1) must be
  72. *          set to .TRUE.; then on exit SELECT(j) is .TRUE. and
  73. *          SELECT(j+1) is .FALSE..
  74. *
  75. *  N       (input) INTEGER
  76. *          The order of the matrix T. N >= 0.
  77. *
  78. *  T       (input) DOUBLE PRECISION array, dimension (LDT,N)
  79. *          The upper quasi-triangular matrix T in Schur canonical form.
  80. *
  81. *  LDT     (input) INTEGER
  82. *          The leading dimension of the array T. LDT >= max(1,N).
  83. *
  84. *  VL      (input/output) DOUBLE PRECISION array, dimension (LDVL,MM)
  85. *          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
  86. *          contain an N-by-N matrix Q (usually the orthogonal matrix Q
  87. *          of Schur vectors returned by DHSEQR).
  88. *          On exit, if SIDE = 'L' or 'B', VL contains:
  89. *          if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
  90. *          if HOWMNY = 'B', the matrix Q*Y;
  91. *          if HOWMNY = 'S', the left eigenvectors of T specified by
  92. *                           SELECT, stored consecutively in the columns
  93. *                           of VL, in the same order as their
  94. *                           eigenvalues.
  95. *          A complex eigenvector corresponding to a complex eigenvalue
  96. *          is stored in two consecutive columns, the first holding the
  97. *          real part, and the second the imaginary part.
  98. *          If SIDE = 'R', VL is not referenced.
  99. *
  100. *  LDVL    (input) INTEGER
  101. *          The leading dimension of the array VL.  LDVL >= max(1,N) if
  102. *          SIDE = 'L' or 'B'; LDVL >= 1 otherwise.
  103. *
  104. *  VR      (input/output) DOUBLE PRECISION array, dimension (LDVR,MM)
  105. *          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
  106. *          contain an N-by-N matrix Q (usually the orthogonal matrix Q
  107. *          of Schur vectors returned by DHSEQR).
  108. *          On exit, if SIDE = 'R' or 'B', VR contains:
  109. *          if HOWMNY = 'A', the matrix X of right eigenvectors of T;
  110. *          if HOWMNY = 'B', the matrix Q*X;
  111. *          if HOWMNY = 'S', the right eigenvectors of T specified by
  112. *                           SELECT, stored consecutively in the columns
  113. *                           of VR, in the same order as their
  114. *                           eigenvalues.
  115. *          A complex eigenvector corresponding to a complex eigenvalue
  116. *          is stored in two consecutive columns, the first holding the
  117. *          real part and the second the imaginary part.
  118. *          If SIDE = 'L', VR is not referenced.
  119. *
  120. *  LDVR    (input) INTEGER
  121. *          The leading dimension of the array VR.  LDVR >= max(1,N) if
  122. *          SIDE = 'R' or 'B'; LDVR >= 1 otherwise.
  123. *
  124. *  MM      (input) INTEGER
  125. *          The number of columns in the arrays VL and/or VR. MM >= M.
  126. *
  127. *  M       (output) INTEGER
  128. *          The number of columns in the arrays VL and/or VR actually
  129. *          used to store the eigenvectors.
  130. *          If HOWMNY = 'A' or 'B', M is set to N.
  131. *          Each selected real eigenvector occupies one column and each
  132. *          selected complex eigenvector occupies two columns.
  133. *
  134. *  WORK    (workspace) DOUBLE PRECISION array, dimension (3*N)
  135. *
  136. *  INFO    (output) INTEGER
  137. *          = 0:  successful exit
  138. *          < 0:  if INFO = -i, the i-th argument had an illegal value
  139. *
  140. *  Further Details
  141. *  ===============
  142. *
  143. *  The algorithm used in this program is basically backward (forward)
  144. *  substitution, with scaling to make the the code robust against
  145. *  possible overflow.
  146. *
  147. *  Each eigenvector is normalized so that the element of largest
  148. *  magnitude has magnitude 1; here the magnitude of a complex number
  149. *  (x,y) is taken to be |x| + |y|.
  150. *
  151. *  =====================================================================
  152. *
  153. *     .. Parameters ..
  154.       DOUBLE PRECISION   ZERO, ONE
  155.       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  156. *     ..
  157. *     .. Local Scalars ..
  158.       LOGICAL            ALLV, BOTHV, LEFTV, OVER, PAIR, RIGHTV, SOMEV
  159.       INTEGER            I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI, N2
  160.       DOUBLE PRECISION   BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE,
  161.      $                   SMIN, SMLNUM, ULP, UNFL, VCRIT, VMAX, WI, WR,
  162.      $                   XNORM
  163. *     ..
  164. *     .. External Functions ..
  165.       LOGICAL            LSAME
  166.       INTEGER            IDAMAX
  167.       DOUBLE PRECISION   DDOT, DLAMCH
  168.       EXTERNAL           LSAME, IDAMAX, DDOT, DLAMCH
  169. *     ..
  170. *     .. External Subroutines ..
  171.       EXTERNAL           DAXPY, DCOPY, DGEMV, DLABAD, DLALN2, DSCAL,
  172.      $                   XERBLA
  173. *     ..
  174. *     .. Intrinsic Functions ..
  175.       INTRINSIC          ABS, MAX, SQRT
  176. *     ..
  177. *     .. Local Arrays ..
  178.       DOUBLE PRECISION   X( 2, 2 )
  179. *     ..
  180. *     .. Executable Statements ..
  181. *
  182. *     Decode and test the input parameters
  183. *
  184.       BOTHV = LSAME( SIDE, 'B' )
  185.       RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
  186.       LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
  187. *
  188.       ALLV = LSAME( HOWMNY, 'A' )
  189.       OVER = LSAME( HOWMNY, 'B' ) .OR. LSAME( HOWMNY, 'O' )
  190.       SOMEV = LSAME( HOWMNY, 'S' )
  191. *
  192.       INFO = 0
  193.       IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
  194.          INFO = -1
  195.       ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
  196.          INFO = -2
  197.       ELSE IF( N.LT.0 ) THEN
  198.          INFO = -4
  199.       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
  200.          INFO = -6
  201.       ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
  202.          INFO = -8
  203.       ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
  204.          INFO = -10
  205.       ELSE
  206. *
  207. *        Set M to the number of columns required to store the selected
  208. *        eigenvectors, standardize the array SELECT if necessary, and
  209. *        test MM.
  210. *
  211.          IF( SOMEV ) THEN
  212.             M = 0
  213.             PAIR = .FALSE.
  214.             DO 10 J = 1, N
  215.                IF( PAIR ) THEN
  216.                   PAIR = .FALSE.
  217.                   SELECT( J ) = .FALSE.
  218.                ELSE
  219.                   IF( J.LT.N ) THEN
  220.                      IF( T( J+1, J ).EQ.ZERO ) THEN
  221.                         IF( SELECT( J ) )
  222.      $                     M = M + 1
  223.                      ELSE
  224.                         PAIR = .TRUE.
  225.                         IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN
  226.                            SELECT( J ) = .TRUE.
  227.                            M = M + 2
  228.                         END IF
  229.                      END IF
  230.                   ELSE
  231.                      IF( SELECT( N ) )
  232.      $                  M = M + 1
  233.                   END IF
  234.                END IF
  235.    10       CONTINUE
  236.          ELSE
  237.             M = N
  238.          END IF
  239. *
  240.          IF( MM.LT.M ) THEN
  241.             INFO = -11
  242.          END IF
  243.       END IF
  244.       IF( INFO.NE.0 ) THEN
  245.          CALL XERBLA( 'DTREVC', -INFO )
  246.          RETURN
  247.       END IF
  248. *
  249. *     Quick return if possible.
  250. *
  251.       IF( N.EQ.0 )
  252.      $   RETURN
  253. *
  254. *     Set the constants to control overflow.
  255. *
  256.       UNFL = DLAMCH( 'Safe minimum' )
  257.       OVFL = ONE / UNFL
  258.       CALL DLABAD( UNFL, OVFL )
  259.       ULP = DLAMCH( 'Precision' )
  260.       SMLNUM = UNFL*( N / ULP )
  261.       BIGNUM = ( ONE-ULP ) / SMLNUM
  262. *
  263. *     Compute 1-norm of each column of strictly upper triangular
  264. *     part of T to control overflow in triangular solver.
  265. *
  266.       WORK( 1 ) = ZERO
  267.       DO 30 J = 2, N
  268.          WORK( J ) = ZERO
  269.          DO 20 I = 1, J - 1
  270.             WORK( J ) = WORK( J ) + ABS( T( I, J ) )
  271.    20    CONTINUE
  272.    30 CONTINUE
  273. *
  274. *     Index IP is used to specify the real or complex eigenvalue:
  275. *       IP = 0, real eigenvalue,
  276. *            1, first of conjugate complex pair: (wr,wi)
  277. *           -1, second of conjugate complex pair: (wr,wi)
  278. *
  279.       N2 = 2*N
  280. *
  281.       IF( RIGHTV ) THEN
  282. *
  283. *        Compute right eigenvectors.
  284. *
  285.          IP = 0
  286.          IS = M
  287.          DO 140 KI = N, 1, -1
  288. *
  289.             IF( IP.EQ.1 )
  290.      $         GO TO 130
  291.             IF( KI.EQ.1 )
  292.      $         GO TO 40
  293.             IF( T( KI, KI-1 ).EQ.ZERO )
  294.      $         GO TO 40
  295.             IP = -1
  296. *
  297.    40       CONTINUE
  298.             IF( SOMEV ) THEN
  299.                IF( IP.EQ.0 ) THEN
  300.                   IF( .NOT.SELECT( KI ) )
  301.      $               GO TO 130
  302.                ELSE
  303.                   IF( .NOT.SELECT( KI-1 ) )
  304.      $               GO TO 130
  305.                END IF
  306.             END IF
  307. *
  308. *           Compute the KI-th eigenvalue (WR,WI).
  309. *
  310.             WR = T( KI, KI )
  311.             WI = ZERO
  312.             IF( IP.NE.0 )
  313.      $         WI = SQRT( ABS( T( KI, KI-1 ) ) )*
  314.      $              SQRT( ABS( T( KI-1, KI ) ) )
  315.             SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
  316. *
  317.             IF( IP.EQ.0 ) THEN
  318. *
  319. *              Real right eigenvector
  320. *
  321.                WORK( KI+N ) = ONE
  322. *
  323. *              Form right-hand side
  324. *
  325.                DO 50 K = 1, KI - 1
  326.                   WORK( K+N ) = -T( K, KI )
  327.    50          CONTINUE
  328. *
  329. *              Solve the upper quasi-triangular system:
  330. *                 (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK.
  331. *
  332.                JNXT = KI - 1
  333.                DO 60 J = KI - 1, 1, -1
  334.                   IF( J.GT.JNXT )
  335.      $               GO TO 60
  336.                   J1 = J
  337.                   J2 = J
  338.                   JNXT = J - 1
  339.                   IF( J.GT.1 ) THEN
  340.                      IF( T( J, J-1 ).NE.ZERO ) THEN
  341.                         J1 = J - 1
  342.                         JNXT = J - 2
  343.                      END IF
  344.                   END IF
  345. *
  346.                   IF( J1.EQ.J2 ) THEN
  347. *
  348. *                    1-by-1 diagonal block
  349. *
  350.                      CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
  351.      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
  352.      $                            ZERO, X, 2, SCALE, XNORM, IERR )
  353. *
  354. *                    Scale X(1,1) to avoid overflow when updating
  355. *                    the right-hand side.
  356. *
  357.                      IF( XNORM.GT.ONE ) THEN
  358.                         IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
  359.                            X( 1, 1 ) = X( 1, 1 ) / XNORM
  360.                            SCALE = SCALE / XNORM
  361.                         END IF
  362.                      END IF
  363. *
  364. *                    Scale if necessary
  365. *
  366.                      IF( SCALE.NE.ONE )
  367.      $                  CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
  368.                      WORK( J+N ) = X( 1, 1 )
  369. *
  370. *                    Update right-hand side
  371. *
  372.                      CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
  373.      $                           WORK( 1+N ), 1 )
  374. *
  375.                   ELSE
  376. *
  377. *                    2-by-2 diagonal block
  378. *
  379.                      CALL DLALN2( .FALSE., 2, 1, SMIN, ONE,
  380.      $                            T( J-1, J-1 ), LDT, ONE, ONE,
  381.      $                            WORK( J-1+N ), N, WR, ZERO, X, 2,
  382.      $                            SCALE, XNORM, IERR )
  383. *
  384. *                    Scale X(1,1) and X(2,1) to avoid overflow when
  385. *                    updating the right-hand side.
  386. *
  387.                      IF( XNORM.GT.ONE ) THEN
  388.                         BETA = MAX( WORK( J-1 ), WORK( J ) )
  389.                         IF( BETA.GT.BIGNUM / XNORM ) THEN
  390.                            X( 1, 1 ) = X( 1, 1 ) / XNORM
  391.                            X( 2, 1 ) = X( 2, 1 ) / XNORM
  392.                            SCALE = SCALE / XNORM
  393.                         END IF
  394.                      END IF
  395. *
  396. *                    Scale if necessary
  397. *
  398.                      IF( SCALE.NE.ONE )
  399.      $                  CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
  400.                      WORK( J-1+N ) = X( 1, 1 )
  401.                      WORK( J+N ) = X( 2, 1 )
  402. *
  403. *                    Update right-hand side
  404. *
  405.                      CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
  406.      $                           WORK( 1+N ), 1 )
  407.                      CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
  408.      $                           WORK( 1+N ), 1 )
  409.                   END IF
  410.    60          CONTINUE
  411. *
  412. *              Copy the vector x or Q*x to VR and normalize.
  413. *
  414.                IF( .NOT.OVER ) THEN
  415.                   CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS ), 1 )
  416. *
  417.                   II = IDAMAX( KI, VR( 1, IS ), 1 )
  418.                   REMAX = ONE / ABS( VR( II, IS ) )
  419.                   CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
  420. *
  421.                   DO 70 K = KI + 1, N
  422.                      VR( K, IS ) = ZERO
  423.    70             CONTINUE
  424.                ELSE
  425.                   IF( KI.GT.1 )
  426.      $               CALL DGEMV( 'N', N, KI-1, ONE, VR, LDVR,
  427.      $                           WORK( 1+N ), 1, WORK( KI+N ),
  428.      $                           VR( 1, KI ), 1 )
  429. *
  430.                   II = IDAMAX( N, VR( 1, KI ), 1 )
  431.                   REMAX = ONE / ABS( VR( II, KI ) )
  432.                   CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
  433.                END IF
  434. *
  435.             ELSE
  436. *
  437. *              Complex right eigenvector.
  438. *
  439. *              Initial solve
  440. *                [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0.
  441. *                [ (T(KI,KI-1)   T(KI,KI)   )               ]
  442. *
  443.                IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN
  444.                   WORK( KI-1+N ) = ONE
  445.                   WORK( KI+N2 ) = WI / T( KI-1, KI )
  446.                ELSE
  447.                   WORK( KI-1+N ) = -WI / T( KI, KI-1 )
  448.                   WORK( KI+N2 ) = ONE
  449.                END IF
  450.                WORK( KI+N ) = ZERO
  451.                WORK( KI-1+N2 ) = ZERO
  452. *
  453. *              Form right-hand side
  454. *
  455.                DO 80 K = 1, KI - 2
  456.                   WORK( K+N ) = -WORK( KI-1+N )*T( K, KI-1 )
  457.                   WORK( K+N2 ) = -WORK( KI+N2 )*T( K, KI )
  458.    80          CONTINUE
  459. *
  460. *              Solve upper quasi-triangular system:
  461. *              (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2)
  462. *
  463.                JNXT = KI - 2
  464.                DO 90 J = KI - 2, 1, -1
  465.                   IF( J.GT.JNXT )
  466.      $               GO TO 90
  467.                   J1 = J
  468.                   J2 = J
  469.                   JNXT = J - 1
  470.                   IF( J.GT.1 ) THEN
  471.                      IF( T( J, J-1 ).NE.ZERO ) THEN
  472.                         J1 = J - 1
  473.                         JNXT = J - 2
  474.                      END IF
  475.                   END IF
  476. *
  477.                   IF( J1.EQ.J2 ) THEN
  478. *
  479. *                    1-by-1 diagonal block
  480. *
  481.                      CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
  482.      $                            LDT, ONE, ONE, WORK( J+N ), N, WR, WI,
  483.      $                            X, 2, SCALE, XNORM, IERR )
  484. *
  485. *                    Scale X(1,1) and X(1,2) to avoid overflow when
  486. *                    updating the right-hand side.
  487. *
  488.                      IF( XNORM.GT.ONE ) THEN
  489.                         IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
  490.                            X( 1, 1 ) = X( 1, 1 ) / XNORM
  491.                            X( 1, 2 ) = X( 1, 2 ) / XNORM
  492.                            SCALE = SCALE / XNORM
  493.                         END IF
  494.                      END IF
  495. *
  496. *                    Scale if necessary
  497. *
  498.                      IF( SCALE.NE.ONE ) THEN
  499.                         CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
  500.                         CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
  501.                      END IF
  502.                      WORK( J+N ) = X( 1, 1 )
  503.                      WORK( J+N2 ) = X( 1, 2 )
  504. *
  505. *                    Update the right-hand side
  506. *
  507.                      CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
  508.      $                           WORK( 1+N ), 1 )
  509.                      CALL DAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1,
  510.      $                           WORK( 1+N2 ), 1 )
  511. *
  512.                   ELSE
  513. *
  514. *                    2-by-2 diagonal block
  515. *
  516.                      CALL DLALN2( .FALSE., 2, 2, SMIN, ONE,
  517.      $                            T( J-1, J-1 ), LDT, ONE, ONE,
  518.      $                            WORK( J-1+N ), N, WR, WI, X, 2, SCALE,
  519.      $                            XNORM, IERR )
  520. *
  521. *                    Scale X to avoid overflow when updating
  522. *                    the right-hand side.
  523. *
  524.                      IF( XNORM.GT.ONE ) THEN
  525.                         BETA = MAX( WORK( J-1 ), WORK( J ) )
  526.                         IF( BETA.GT.BIGNUM / XNORM ) THEN
  527.                            REC = ONE / XNORM
  528.                            X( 1, 1 ) = X( 1, 1 )*REC
  529.                            X( 1, 2 ) = X( 1, 2 )*REC
  530.                            X( 2, 1 ) = X( 2, 1 )*REC
  531.                            X( 2, 2 ) = X( 2, 2 )*REC
  532.                            SCALE = SCALE*REC
  533.                         END IF
  534.                      END IF
  535. *
  536. *                    Scale if necessary
  537. *
  538.                      IF( SCALE.NE.ONE ) THEN
  539.                         CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
  540.                         CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
  541.                      END IF
  542.                      WORK( J-1+N ) = X( 1, 1 )
  543.                      WORK( J+N ) = X( 2, 1 )
  544.                      WORK( J-1+N2 ) = X( 1, 2 )
  545.                      WORK( J+N2 ) = X( 2, 2 )
  546. *
  547. *                    Update the right-hand side
  548. *
  549.                      CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
  550.      $                           WORK( 1+N ), 1 )
  551.                      CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
  552.      $                           WORK( 1+N ), 1 )
  553.                      CALL DAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1,
  554.      $                           WORK( 1+N2 ), 1 )
  555.                      CALL DAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1,
  556.      $                           WORK( 1+N2 ), 1 )
  557.                   END IF
  558.    90          CONTINUE
  559. *
  560. *              Copy the vector x or Q*x to VR and normalize.
  561. *
  562.                IF( .NOT.OVER ) THEN
  563.                   CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS-1 ), 1 )
  564.                   CALL DCOPY( KI, WORK( 1+N2 ), 1, VR( 1, IS ), 1 )
  565. *
  566.                   EMAX = ZERO
  567.                   DO 100 K = 1, KI
  568.                      EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+
  569.      $                      ABS( VR( K, IS ) ) )
  570.   100             CONTINUE
  571. *
  572.                   REMAX = ONE / EMAX
  573.                   CALL DSCAL( KI, REMAX, VR( 1, IS-1 ), 1 )
  574.                   CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
  575. *
  576.                   DO 110 K = KI + 1, N
  577.                      VR( K, IS-1 ) = ZERO
  578.                      VR( K, IS ) = ZERO
  579.   110             CONTINUE
  580. *
  581.                ELSE
  582. *
  583.                   IF( KI.GT.2 ) THEN
  584.                      CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
  585.      $                           WORK( 1+N ), 1, WORK( KI-1+N ),
  586.      $                           VR( 1, KI-1 ), 1 )
  587.                      CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
  588.      $                           WORK( 1+N2 ), 1, WORK( KI+N2 ),
  589.      $                           VR( 1, KI ), 1 )
  590.                   ELSE
  591.                      CALL DSCAL( N, WORK( KI-1+N ), VR( 1, KI-1 ), 1 )
  592.                      CALL DSCAL( N, WORK( KI+N2 ), VR( 1, KI ), 1 )
  593.                   END IF
  594. *
  595.                   EMAX = ZERO
  596.                   DO 120 K = 1, N
  597.                      EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+
  598.      $                      ABS( VR( K, KI ) ) )
  599.   120             CONTINUE
  600.                   REMAX = ONE / EMAX
  601.                   CALL DSCAL( N, REMAX, VR( 1, KI-1 ), 1 )
  602.                   CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
  603.                END IF
  604.             END IF
  605. *
  606.             IS = IS - 1
  607.             IF( IP.NE.0 )
  608.      $         IS = IS - 1
  609.   130       CONTINUE
  610.             IF( IP.EQ.1 )
  611.      $         IP = 0
  612.             IF( IP.EQ.-1 )
  613.      $         IP = 1
  614.   140    CONTINUE
  615.       END IF
  616. *
  617.       IF( LEFTV ) THEN
  618. *
  619. *        Compute left eigenvectors.
  620. *
  621.          IP = 0
  622.          IS = 1
  623.          DO 260 KI = 1, N
  624. *
  625.             IF( IP.EQ.-1 )
  626.      $         GO TO 250
  627.             IF( KI.EQ.N )
  628.      $         GO TO 150
  629.             IF( T( KI+1, KI ).EQ.ZERO )
  630.      $         GO TO 150
  631.             IP = 1
  632. *
  633.   150       CONTINUE
  634.             IF( SOMEV ) THEN
  635.                IF( .NOT.SELECT( KI ) )
  636.      $            GO TO 250
  637.             END IF
  638. *
  639. *           Compute the KI-th eigenvalue (WR,WI).
  640. *
  641.             WR = T( KI, KI )
  642.             WI = ZERO
  643.             IF( IP.NE.0 )
  644.      $         WI = SQRT( ABS( T( KI, KI+1 ) ) )*
  645.      $              SQRT( ABS( T( KI+1, KI ) ) )
  646.             SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
  647. *
  648.             IF( IP.EQ.0 ) THEN
  649. *
  650. *              Real left eigenvector.
  651. *
  652.                WORK( KI+N ) = ONE
  653. *
  654. *              Form right-hand side
  655. *
  656.                DO 160 K = KI + 1, N
  657.                   WORK( K+N ) = -T( KI, K )
  658.   160          CONTINUE
  659. *
  660. *              Solve the quasi-triangular system:
  661. *                 (T(KI+1:N,KI+1:N) - WR)'*X = SCALE*WORK
  662. *
  663.                VMAX = ONE
  664.                VCRIT = BIGNUM
  665. *
  666.                JNXT = KI + 1
  667.                DO 170 J = KI + 1, N
  668.                   IF( J.LT.JNXT )
  669.      $               GO TO 170
  670.                   J1 = J
  671.                   J2 = J
  672.                   JNXT = J + 1
  673.                   IF( J.LT.N ) THEN
  674.                      IF( T( J+1, J ).NE.ZERO ) THEN
  675.                         J2 = J + 1
  676.                         JNXT = J + 2
  677.                      END IF
  678.                   END IF
  679. *
  680.                   IF( J1.EQ.J2 ) THEN
  681. *
  682. *                    1-by-1 diagonal block
  683. *
  684. *                    Scale if necessary to avoid overflow when forming
  685. *                    the right-hand side.
  686. *
  687.                      IF( WORK( J ).GT.VCRIT ) THEN
  688.                         REC = ONE / VMAX
  689.                         CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
  690.                         VMAX = ONE
  691.                         VCRIT = BIGNUM
  692.                      END IF
  693. *
  694.                      WORK( J+N ) = WORK( J+N ) -
  695.      $                             DDOT( J-KI-1, T( KI+1, J ), 1,
  696.      $                             WORK( KI+1+N ), 1 )
  697. *
  698. *                    Solve (T(J,J)-WR)'*X = WORK
  699. *
  700.                      CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
  701.      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
  702.      $                            ZERO, X, 2, SCALE, XNORM, IERR )
  703. *
  704. *                    Scale if necessary
  705. *
  706.                      IF( SCALE.NE.ONE )
  707.      $                  CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
  708.                      WORK( J+N ) = X( 1, 1 )
  709.                      VMAX = MAX( ABS( WORK( J+N ) ), VMAX )
  710.                      VCRIT = BIGNUM / VMAX
  711. *
  712.                   ELSE
  713. *
  714. *                    2-by-2 diagonal block
  715. *
  716. *                    Scale if necessary to avoid overflow when forming
  717. *                    the right-hand side.
  718. *
  719.                      BETA = MAX( WORK( J ), WORK( J+1 ) )
  720.                      IF( BETA.GT.VCRIT ) THEN
  721.                         REC = ONE / VMAX
  722.                         CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
  723.                         VMAX = ONE
  724.                         VCRIT = BIGNUM
  725.                      END IF
  726. *
  727.                      WORK( J+N ) = WORK( J+N ) -
  728.      $                             DDOT( J-KI-1, T( KI+1, J ), 1,
  729.      $                             WORK( KI+1+N ), 1 )
  730. *
  731.                      WORK( J+1+N ) = WORK( J+1+N ) -
  732.      $                               DDOT( J-KI-1, T( KI+1, J+1 ), 1,
  733.      $                               WORK( KI+1+N ), 1 )
  734. *
  735. *                    Solve
  736. *                      [T(J,J)-WR   T(J,J+1)     ]'* X = SCALE*( WORK1 )
  737. *                      [T(J+1,J)    T(J+1,J+1)-WR]             ( WORK2 )
  738. *
  739.                      CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ),
  740.      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
  741.      $                            ZERO, X, 2, SCALE, XNORM, IERR )
  742. *
  743. *                    Scale if necessary
  744. *
  745.                      IF( SCALE.NE.ONE )
  746.      $                  CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
  747.                      WORK( J+N ) = X( 1, 1 )
  748.                      WORK( J+1+N ) = X( 2, 1 )
  749. *
  750.                      VMAX = MAX( ABS( WORK( J+N ) ),
  751.      $                      ABS( WORK( J+1+N ) ), VMAX )
  752.                      VCRIT = BIGNUM / VMAX
  753. *
  754.                   END IF
  755.   170          CONTINUE
  756. *
  757. *              Copy the vector x or Q*x to VL and normalize.
  758. *
  759.                IF( .NOT.OVER ) THEN
  760.                   CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
  761. *
  762.                   II = IDAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
  763.                   REMAX = ONE / ABS( VL( II, IS ) )
  764.                   CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
  765. *
  766.                   DO 180 K = 1, KI - 1
  767.                      VL( K, IS ) = ZERO
  768.   180             CONTINUE
  769. *
  770.                ELSE
  771. *
  772.                   IF( KI.LT.N )
  773.      $               CALL DGEMV( 'N', N, N-KI, ONE, VL( 1, KI+1 ), LDVL,
  774.      $                           WORK( KI+1+N ), 1, WORK( KI+N ),
  775.      $                           VL( 1, KI ), 1 )
  776. *
  777.                   II = IDAMAX( N, VL( 1, KI ), 1 )
  778.                   REMAX = ONE / ABS( VL( II, KI ) )
  779.                   CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
  780. *
  781.                END IF
  782. *
  783.             ELSE
  784. *
  785. *              Complex left eigenvector.
  786. *
  787. *               Initial solve:
  788. *                 ((T(KI,KI)    T(KI,KI+1) )' - (WR - I* WI))*X = 0.
  789. *                 ((T(KI+1,KI) T(KI+1,KI+1))                )
  790. *
  791.                IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN
  792.                   WORK( KI+N ) = WI / T( KI, KI+1 )
  793.                   WORK( KI+1+N2 ) = ONE
  794.                ELSE
  795.                   WORK( KI+N ) = ONE
  796.                   WORK( KI+1+N2 ) = -WI / T( KI+1, KI )
  797.                END IF
  798.                WORK( KI+1+N ) = ZERO
  799.                WORK( KI+N2 ) = ZERO
  800. *
  801. *              Form right-hand side
  802. *
  803.                DO 190 K = KI + 2, N
  804.                   WORK( K+N ) = -WORK( KI+N )*T( KI, K )
  805.                   WORK( K+N2 ) = -WORK( KI+1+N2 )*T( KI+1, K )
  806.   190          CONTINUE
  807. *
  808. *              Solve complex quasi-triangular system:
  809. *              ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2
  810. *
  811.                VMAX = ONE
  812.                VCRIT = BIGNUM
  813. *
  814.                JNXT = KI + 2
  815.                DO 200 J = KI + 2, N
  816.                   IF( J.LT.JNXT )
  817.      $               GO TO 200
  818.                   J1 = J
  819.                   J2 = J
  820.                   JNXT = J + 1
  821.                   IF( J.LT.N ) THEN
  822.                      IF( T( J+1, J ).NE.ZERO ) THEN
  823.                         J2 = J + 1
  824.                         JNXT = J + 2
  825.                      END IF
  826.                   END IF
  827. *
  828.                   IF( J1.EQ.J2 ) THEN
  829. *
  830. *                    1-by-1 diagonal block
  831. *
  832. *                    Scale if necessary to avoid overflow when
  833. *                    forming the right-hand side elements.
  834. *
  835.                      IF( WORK( J ).GT.VCRIT ) THEN
  836.                         REC = ONE / VMAX
  837.                         CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
  838.                         CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
  839.                         VMAX = ONE
  840.                         VCRIT = BIGNUM
  841.                      END IF
  842. *
  843.                      WORK( J+N ) = WORK( J+N ) -
  844.      $                             DDOT( J-KI-2, T( KI+2, J ), 1,
  845.      $                             WORK( KI+2+N ), 1 )
  846.                      WORK( J+N2 ) = WORK( J+N2 ) -
  847.      $                              DDOT( J-KI-2, T( KI+2, J ), 1,
  848.      $                              WORK( KI+2+N2 ), 1 )
  849. *
  850. *                    Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2
  851. *
  852.                      CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
  853.      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
  854.      $                            -WI, X, 2, SCALE, XNORM, IERR )
  855. *
  856. *                    Scale if necessary
  857. *
  858.                      IF( SCALE.NE.ONE ) THEN
  859.                         CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
  860.                         CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
  861.                      END IF
  862.                      WORK( J+N ) = X( 1, 1 )
  863.                      WORK( J+N2 ) = X( 1, 2 )
  864.                      VMAX = MAX( ABS( WORK( J+N ) ),
  865.      $                      ABS( WORK( J+N2 ) ), VMAX )
  866.                      VCRIT = BIGNUM / VMAX
  867. *
  868.                   ELSE
  869. *
  870. *                    2-by-2 diagonal block
  871. *
  872. *                    Scale if necessary to avoid overflow when forming
  873. *                    the right-hand side elements.
  874. *
  875.                      BETA = MAX( WORK( J ), WORK( J+1 ) )
  876.                      IF( BETA.GT.VCRIT ) THEN
  877.                         REC = ONE / VMAX
  878.                         CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
  879.                         CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
  880.                         VMAX = ONE
  881.                         VCRIT = BIGNUM
  882.                      END IF
  883. *
  884.                      WORK( J+N ) = WORK( J+N ) -
  885.      $                             DDOT( J-KI-2, T( KI+2, J ), 1,
  886.      $                             WORK( KI+2+N ), 1 )
  887. *
  888.                      WORK( J+N2 ) = WORK( J+N2 ) -
  889.      $                              DDOT( J-KI-2, T( KI+2, J ), 1,
  890.      $                              WORK( KI+2+N2 ), 1 )
  891. *
  892.                      WORK( J+1+N ) = WORK( J+1+N ) -
  893.      $                               DDOT( J-KI-2, T( KI+2, J+1 ), 1,
  894.      $                               WORK( KI+2+N ), 1 )
  895. *
  896.                      WORK( J+1+N2 ) = WORK( J+1+N2 ) -
  897.      $                                DDOT( J-KI-2, T( KI+2, J+1 ), 1,
  898.      $                                WORK( KI+2+N2 ), 1 )
  899. *
  900. *                    Solve 2-by-2 complex linear equation
  901. *                      ([T(j,j)   T(j,j+1)  ]'-(wr-i*wi)*I)*X = SCALE*B
  902. *                      ([T(j+1,j) T(j+1,j+1)]             )
  903. *
  904.                      CALL DLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ),
  905.      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
  906.      $                            -WI, X, 2, SCALE, XNORM, IERR )
  907. *
  908. *                    Scale if necessary
  909. *
  910.                      IF( SCALE.NE.ONE ) THEN
  911.                         CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
  912.                         CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
  913.                      END IF
  914.                      WORK( J+N ) = X( 1, 1 )
  915.                      WORK( J+N2 ) = X( 1, 2 )
  916.                      WORK( J+1+N ) = X( 2, 1 )
  917.                      WORK( J+1+N2 ) = X( 2, 2 )
  918.                      VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ),
  919.      $                      ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ), VMAX )
  920.                      VCRIT = BIGNUM / VMAX
  921. *
  922.                   END IF
  923.   200          CONTINUE
  924. *
  925. *              Copy the vector x or Q*x to VL and normalize.
  926. *
  927.   210          CONTINUE
  928.                IF( .NOT.OVER ) THEN
  929.                   CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
  930.                   CALL DCOPY( N-KI+1, WORK( KI+N2 ), 1, VL( KI, IS+1 ),
  931.      $                        1 )
  932. *
  933.                   EMAX = ZERO
  934.                   DO 220 K = KI, N
  935.                      EMAX = MAX( EMAX, ABS( VL( K, IS ) )+
  936.      $                      ABS( VL( K, IS+1 ) ) )
  937.   220             CONTINUE
  938.                   REMAX = ONE / EMAX
  939.                   CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
  940.                   CALL DSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 )
  941. *
  942.                   DO 230 K = 1, KI - 1
  943.                      VL( K, IS ) = ZERO
  944.                      VL( K, IS+1 ) = ZERO
  945.   230             CONTINUE
  946.                ELSE
  947.                   IF( KI.LT.N-1 ) THEN
  948.                      CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
  949.      $                           LDVL, WORK( KI+2+N ), 1, WORK( KI+N ),
  950.      $                           VL( 1, KI ), 1 )
  951.                      CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
  952.      $                           LDVL, WORK( KI+2+N2 ), 1,
  953.      $                           WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
  954.                   ELSE
  955.                      CALL DSCAL( N, WORK( KI+N ), VL( 1, KI ), 1 )
  956.                      CALL DSCAL( N, WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
  957.                   END IF
  958. *
  959.                   EMAX = ZERO
  960.                   DO 240 K = 1, N
  961.                      EMAX = MAX( EMAX, ABS( VL( K, KI ) )+
  962.      $                      ABS( VL( K, KI+1 ) ) )
  963.   240             CONTINUE
  964.                   REMAX = ONE / EMAX
  965.                   CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
  966.                   CALL DSCAL( N, REMAX, VL( 1, KI+1 ), 1 )
  967. *
  968.                END IF
  969. *
  970.             END IF
  971. *
  972.             IS = IS + 1
  973.             IF( IP.NE.0 )
  974.      $         IS = IS + 1
  975.   250       CONTINUE
  976.             IF( IP.EQ.-1 )
  977.      $         IP = 0
  978.             IF( IP.EQ.1 )
  979.      $         IP = -1
  980. *
  981.   260    CONTINUE
  982. *
  983.       END IF
  984. *
  985.       RETURN
  986. *
  987. *     End of DTREVC
  988. *
  989.       END
  990.